var Yt Ct Bt Lpt Kt Nt Wt Rt Qt Zt Gt Ot Mut;

parameters beta varphi delta psi F nu alpha theta phi;


load xpar;
 
beta       = xpar(1);
varphi     = xpar(2);
delta      = xpar(3);
psi        = xpar(4);
F          = xpar(5);
nu         = xpar(6);
alpha      = xpar(7);
theta      = xpar(8);
phi        = xpar(9);

external_function(name = interpv,  nargs = 2, first_deriv_provided, second_deriv_provided); 
 

model;

1          = beta*Ct/Ct(+1)*(Rt(+1) + 1 - delta);
F*Wt       = beta*Ct/Ct(+1)*Qt(+1);
Qt         = Yt/Nt(-1)*Gt + beta*(1 - varphi)*Ct/Ct(+1)*Qt(+1);
Wt         = psi*Ct*(Lpt + F*(Nt - (1 - varphi)*Nt(-1)))^nu; 
Ot         = (phi*((Rt/alpha)^alpha*(Wt/(1 - alpha))^(1 - alpha))^(1 - theta) + (1 - phi))^(1/(1 - theta)); 
Yt         = Ct + Kt - (1 - delta)*Kt(-1) + Bt; 
Wt         = (1 - alpha)*phi/Mut*((Rt/alpha)^alpha*(Wt/(1 - alpha))^(1 - alpha)/Ot)^(1 - theta)*Yt/Lpt; 
Rt         =       alpha*phi/Mut*((Rt/alpha)^alpha*(Wt/(1 - alpha))^(1 - alpha)/Ot)^(1 - theta)*Yt/Kt(-1);
1          = (1 - phi)/Mut*(1/Ot)^(1 - theta)*Yt/Bt;

Zt         = interpv(Nt(-1), 1);
Gt         = interpv(Nt(-1), 2);
Mut        = interpv(Nt(-1), 3);

Ot         = Zt/Mut;


end;

load x1_ss;
load x2_ss;


initval;

Kt   = x1_ss(1); 
Nt   = x1_ss(2);

end;


endval;


Yt   = x2_ss(1);
Ct   = x2_ss(2);
Bt   = x2_ss(3);
Lpt  = x2_ss(4);
Kt   = x2_ss(5);
Nt   = x2_ss(6);
Wt   = x2_ss(7);
Rt   = x2_ss(8);
Qt   = x2_ss(9);
Zt   = x2_ss(10);
Gt   = x2_ss(11);
Ot   = x2_ss(12);
Mut  = x2_ss(13);  


end;



resid;
steady;

perfect_foresight_setup(periods = 500);
perfect_foresight_solver(tolf = 1e-10, tolx = 1e-10);
